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ABSTRACT 

Quasi-spherical supersonic (Bondi-type) accretion to a star with a dipole magnetic field is investigated 
using resistive magnetohydrodynamic simulations. A systematic study is made of accretion to a non- 
rotating star, while sample results for a rotating star are also presented. We find that an approximately 
spherical shock wave forms around the dipole with an essential part of the star's initial magnetic flux 
compressed inside the shock wave. A new stationary subsonic accretion flow is established inside the 
shock wave with a steady rate of accretion to the star smaller than the Bondi accretion rate Mb- Matter 
accumulates between the star and the shock wave with the result that the shock wave expands. Accretion 
to the dipole is almost spherically symmetric at radii larger than 2Ra, where Ra is the Alfven radius, 
but it is strongly anisotropic at distances comparable to the Alfven radius and smaller. At these small 
distances matter flows along the magnetic field lines and accretes to the poles of the star along polar 
columns. The accretion Sow becomes supersonic in the region of the polar columns. In a test case with 
an unmagnetized star, we observed spherically-symmetric stationary Bondi accretion without a shock 
wave. The accretion rate to the dipole Mdip is found to depend on j3 oc Ms/fJ. 2 , where fj, is the star's 
magnetic moment, and r\ m the magnetic diffusivity. Specifically, Mdi P oc /3 ' 5 and Mdi P oc ?7„ 38 - The 
equatorial Alfven radius is found to depend on f3 as Ra °c /3 -0 ' 3 which is close to theoretical dependence 
oc /3" 2 / 7 . There is a weak dependence on magnetic diffusivity, Ra oc 77^ . 

Simulations of accretion to a rotating star with an aligned dipole magnetic field show that for slow 
rotation the accretion flow is similar to that in non-rotating case with somewhat smaller values of Mdi p . 
In the case of fast rotation the structure of the subsonic accretion flow is fundamentally different and 
includes a region of "propeller" outflow. The methods and results described here are of general interest 
and can be applied to systems where matter accretes with low angular momentum. 

Subject headings: accretion, dipole — plasmas — magnetic fields — stars: magnetic fields — X-rays: 
stars 



1. INTRODUCTION 

Accretion of matter to a rotating star with a dipole mag- 
netic field is a complex and still unsolved problem in as- 
trophysics. The simplest limit is that of accretion to a star 
with an aligned dipole magnetic field. Although in many 
cases accretion occurs through a disk, in other cases, where 
accreting matter has small angular momentum the accre- 
tion flow is quasi-spherical at large distances from the star. 
Examples include some types of wind fed pulsars (see re- 
view by Nagase 1989). Also, quasi-spherical accretion may 
occur to an isolated star if its velocity through the inter- 
stellar medium is small compared with the sound speed. 
Advection dominated accretion ( Paczyhski & Bisnovatyi- 



Kogan 1981; Narayan & Yi 1995) is also expected to be 
quasi-spherical. 

A general analytic solution for spherical accretion to a 
non-magnetized star was obtained by Bondi (1952). His 
results have also been confirmed now by numerical three- 
dimensional (3D) hydrodynamic simulations by Ruffert 
(1994). The theory and simulations show that matter ac- 
cretes steadily to the gravitating center without formation 
of shocks. Accretion of matter with low angular momen- 
tum to non-magnetized center was investigated recently by 
Bisnovatyi-Kogan & Pogorelov (1997). Less attention has 
been given to quasi-spherical accretion to a magnetized 
star. Disk accretion to a rotating star with an aligned 
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dipole magnetic field has been investigated in a number 
of papers both analytically (Pringle & Rees 1972; Ghosh 
& Lamb 1978; Wang 1979; Shu et al. 1988; Lovelace, 
Romanova, & Bisnovatyi-Kogan 1995, 1998; Li & Wick- 
ramasinghe 1997) and by numerical simulations (Hayashi, 
Shibata, & Matsumoto 1996; Goodson, Winglee, & Bohm 
1997; Miller & Stone 1997). 

Investigation of quasi-spherical accretion to a rotating 
star with dipole field is important because it is a relatively 
simple limit where different aspects of accretion to a dipole 
can be observed and clarified. The general nature of quasi- 
spherical accretion was proposed earlier (Davidson & Os- 
triker 1973; Lamb, Pethick, & Pines 1973; Arons & Lea 
1976; Lipunov 1992, and references therein), but the the- 
oretical ideas have not been tested by MHD simulations. 
The questions of interest include the global nature of the 
accretion flow, the location and the shape of the Alfven 
surface, and the flow structure, in particular, the depar- 
tures of the flow from spherical inflow to highly anisotropic 
polar column accretion inside the dipole's magnetosphere. 
Also, it is of interest to verify dependence of the Alfven 
radius on the accretion rate, the star's magnetic moment 
and rotation rate, and the magnetic diffusivity (considered 
by Lovelace et al. 1995 for the case of disk accretion). 

This paper investigates spherical accretion to a rotating 
star with an aligned dipole magnetic field by axisymmet- 
ric, time-dependent, resistive MHD simulations. Section 
2 describes the model, the equations, the boundary and 
initial conditions, and the numerical methods used. Sec- 
tion 3 discusses the results of simulations for non-rotating 
and rotating central object. A numerical astrophysical ex- 
ample is given in §3.5. Section 4 gives the conclusions of 
this work. 

2. MODEL 

Here, we describe the approach we have taken in ax- 
isymmetric MHD simulations of accretion to a rotating 
star with an aligned dipole magnetic field. We present 
the mathematical model, including the complete system 
of resistive MHD equations, the method used to establish 
the star's intrinsic dipole magnetic field, the initial and 
boundary conditions, and a description of the numerical 
method used to solve the MHD equations. 

2.1. System of Equations 

We consider the equation system for resistive MHD 
(Landau & Lifshitz 1960), 
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All variables have their usual meanings. The equation 
of state is considered to be that for an ideal gas, p = 
(7 — 1) pe, with 7 the usual specific heat ratio. The equa- 
tions incorporate Ohm's law J = cr(E + v x H/c), where 



a is the electrical conductivity. In equation (2) the gravi- 
tational force F 9 (R) = —GMpH/R 3 , is due to the central 
star, where R is the radius vector, and M is the star's 
mass. 

We use an inertial cylindrical coordinate system (r, (f>, z). 
The z-axis is parallel to the star's rotation axis and dipole 
magnetic moment fi. The coordinate system origin coin- 
cides with the star's center and dipole's center. Axisym- 
metry is assumed, d/d(f> — 0. Further, symmetry about 
the z = plane is assumed. Thus calculations may be 
performed on one-quarter of the r — z plane so that the 
"computational region" is < R < R ma x, < z < Z max . 
A totally absorbing sphere, an "accretor" was placed close 
around the origin. The radius of the accretor was chosen 
to be small, r accr < R max - 

In order to guarantee that V • H = holds for all time 
in the numerical simulations, we use the vector potential 
A for the magnetic field, H = V x A, instead of magnetic 
field H itself. For axisymmetric conditions equations (1) 
- (4) can be written in terms of the toroidal vector poten- 
tial (or the flux function * = rA^) and of the toroidal 
magnetic field H^: 
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Here, we have introduced the magnetic diffusivity r\ m = 
c 2 /(A-ko-) — const, and = pv^r, which is the angular 
momentum density. The poloidal components of the mag- 
netic field are H r = —dA^/dz and H z = (l/r)d(rA^)/dr. 

2.2. Method of Establishing Star's Dipole Field 

The intrinsic magnetic field of the star is generated by 
current-density J flowing inside it. In the absence of 
plasma currents outside of the star, the vector potential at 
a point R is A (R) = (1/c) J d 3 x' J(R)/|R-R'|. At large 
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Fig. 1. — The figure shows the geometry of the model. A rotating star of mass M and radius R is replaced by in- 
finitcsimely thin, rotating, "superconducting" disk (0 < r < Rd,z = 0) which is a part of the boundary condition. A 
fixed azimuthal current flows in this disk and generates the intrinsic dipole-likc magnetic field of the star. An absorbing 
sphere, or "an accretor" with radius R aC c = O.5R4 is located at the center. We assume axisymmetry and reflection sym- 
metry about equatorial plane. Thus calculations need to be performed only in one quarter of the (r, z) plane, in the box 
(0 < r < Rmax, < z < Z max ). A spherically symmetric supersonic infow of matter occurs through the outer boundaries 
of the computational region. The inflow rate is Bondi rate Mb as discussed in the text. 



distances from the star, the vector potential can be approx- 
imated as A = fi x R/i? 3 , where fi = (1/c) J d 3 x R x J 
is the intrinsic magnetic moment of the star. The cor- 
responding magnetic field is H = [3r Qu • r) — r 2 |/x|]/r 5 , 
which is a "pure" dipole field. 

In order to establish an intrinsic stellar dipole field in 
our simulations we introduce an "external" surface cur- 
rent flowing on a finite part of the equatorial plane, that 
is, in a disk in the region < Rd << Rmax- This current 
models the current flowing inside the star. There are no 
additional external currents in our model. The presence 
of this "current disk" creates a dipole-type intrinsic mag- 
netic field in our computational box. The nature of this 
field is shown in Figure [l]. 

We choose the azimuthal current-density of the "current 
disk" to be 



il:J,Ar. : ) J n ( — j 



r 



(12) 



for < r < Rd <C Rmax and z = 0, where j\ and j'2 are 
constants. The magnetic moment of this current is 



Rd 

A* = J dr'r' 2 j<f, [r') 
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For the current distribution (12) with j\ = 3 and j% = 1 
(used subsequently in our simulations), 



fi = e 
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The vector potential corresponding to the azimuthal cur- 
rent density (12) at R = (r, z) is 



A4, ( r ' z) = - dr — \ — 
c J k V r 
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where k 2 = Arr' / (r + r') 2 + z 



, and K, E are the full 



elliptic integrals of the first and the second type, 
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We use equations (12) and (15) to numerically determine 
Ac/, in the computational region including the surface of 
the "current disk" . 

The vacuum magnetic field of the "current disk" is given 
by H = V x (A^cj)) with A^ given by (15). For example, 
the magnetic field at the center of the disk for j\ = 3 and 
.72 = 1 is 

ttJq 7fi 



H(0,0)=e, 
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The magnetic field is found to be close to that of a point 
dipole for R > 1.5Rd- The initial (vacuum) magnetic field 
is shown at Figure [j]. 

2.3. Boundary and Initial Conditions 

Here, we consider the boundary conditions on the four 
sides of our computational region < r < R ma x , < z < 
Zmax (see Figure |l]). We first consider the conditions on 
the bottom boundary (0 < r < Rmax-, 2 = 0). The region 
of the above mentioned "current disk" (0 < r < Rdisk < 
Rmax, z = 0) we treat as perfectly conducting or in effect 
"superconducting" . We consider the general case where 
this disk, which represents the star, rotates rigidly with 
angular rate ui about the z— axis. Consequently, the elec- 
tric field in the comoving frame of the disk is zero. In the 
laboratory or non-rotating reference frame, the tangential 
components of the electric field at (0 < r < Rdisk, z = 0) 
are 

1 



E r (r,0) 
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These relations hold for all time in our simulations. 

Owing to the assumed axisymmctry, E 
-(l/c)dA/dt- V$ = -(l/c)dA/dt so that dA+/dt 
—cE^. Consequently, equation (17) gives 



dA^r,0) 
dt 



= (0 < r < R d ) 



(18) 



We also have 

7 disk _ Vm dH,/, 
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c dz 

or, using equation (17), 
(v,p - ur)H z (r, 0) = r/ r 



dz 



{v<pH z - v z H,f,) , 



(0 < r < R d ) . (19) 



For simplicity we consider that the radial current den- 
sity is zero at the disk for < r < Rdisk, Jr = 
— [c/iAir^dH^/dz — 0. From this we obtain the condition 
v<p = lot on the current disk part of the z = boundary. 
In effect we have a no slip condition on the current disk. 
The full set of boundary conditions on the current disk 
(0 < r < R dl z = 0) are 



v z (r,0) = Q, v r (r,0) = 0, u*(r,0)=wr, 



dA^ (r,0) 
dt 







dz 







(20) 



The condition dA^ (r, 0)/ dt = implies that the vector 
potential at the surface of the current disk is independent 
of time. The potential A^ (r, 0) on this surface is obtained 
at the beginning of the simulation using equation (15), and 
it is fixed during the simulation. 

The region of the equatorial plane outside of the cur- 
rent disk {Rd < r < R ma x , z = 0) is treated as a symmetry 
plane. Thus in this region the boundary conditions are 



dA^(r,0) 
dz 







H 4 ,(r,0) = 0, v z (r,0)=0. (21) 



Owing to the assumed axisymmetry, the boundary con- 
ditions on the z-axis (r = , < z < Z max ) are 



v r (0,z) = , Vj,(0,z) = , 
^(0,z) = 0, H <t> (0,z) = 0. 



(22) 



Next, we consider the conditions at the outer bound- 
aries. For these boundaries we assume spherically sym- 
metrical inflow with physical values given by the classical 
Bondi (1952) solution (see also Holzer & Axford 1970). 
The accretion rate 



M = 4ttA 



(f)" 



(23) 



in the Bondi solution is defined by the density p^ and 
the sound speed Coo at infinity and by the mass of the 
central object M. The characteristic length-scale of the 
problem is the Bondi radius Rb = GM/c^. The type 
of solution is defined by the value of the dimensionless 
parameter A. The maximum possible value of A, de- 
noted A c , (for example, A c = 0.625 for 7 = 7/5), cor- 
responds to a flow which is subsonic at infinity and su- 
personic within the sonic point at a radius r s — [(5 — 



3^)/4]Rb. We used the maximum A Bondi solution for 
given parameters at infinity for establishing the conditions 
at the inflow boundaries (r — R m ax, < z < Z max ) and 
(0 < r < Rmax, z — Z max ). The maximum accretion rate 
referred to as Mb correponds to the maximum A. Because 
the inflow is supersonic, all gas dynamical variables can be 
fixed at the outer boundary. 

At the outer inflow boundaries, the magnetic field is 
assumed to vanish, = 0, A^ = 0. This is reasonable 
for situations where the inflowing plasma is unmagnetized. 
Such an inflow will tend to "carry" inward the dipole field 
of the star which is relatively weak at the outer bound- 
aries. However, for comparison we tested different condi- 
tions for non-zero A^ on the outer boundaries, for example, 
dA^/dn = 0. Our numerical results show no dependence 
of the flow on the outer conditions on A^. 

An absorbing sphere or an "accretor" is located close 
around the origin (r, z) = 0. It prevents matter accumu- 
lation in the region close to the gravitating center. The 
radius of the accretor is taken to be R aC cr — 0.5R d (see 
Figure [j]). After each time step, the matter pressure inside 
this sphere is set to a small value in comparison with the 
pressure in the immediately surrounding region. The typ- 
ical pressure contrast was 1/1000. This allows the matter 
just outside the "accretor" to expand freely inward into the 
low pressure region R < R aC cr- Thus, super slowmagne- 
tosonic inflow into the accretor is realized during the sim- 
ulations. In contrast to pure hydrodynamic simulations 
(for example, Ruffert 1994), the density is not small in- 
side the "accretor" . Close to the origin, the magnetic field 
has its highest value. A low density in this region gives a 
high Alfven speed and therefore a very small time step as 
follows from the Courant-Friedrichs-Levy condition. For 
this reason the density inside the "accretor" was set equal 
to a fraction of the exterior density while the temperature 
was set to a small fraction of the exterior temperature in 
order to provide the mentioned pressure contrast with the 
matter just outside the accretor. 

At t = 0, magnetic field is obtained from the vacuum 
vector potential A^ of equation (15). For all runs, with 
non-rotating and rotating central objects, the azimuthal 
magnetic field is zero, H^r.t = 0) = 0. Density, pressure, 
and velocity fields were taken from the Bondi solution with 
maximum possible accretion rate. Inside a sphere with ra- 
dius R = 3R d , the velocity at t = was set to zero. 



2.4. Dimensionless Parameters and Variables 

It is helpful to put equations (5)— (11) into dimensionless 
form. For this we use the values poo, Hq, and R d for the 
density, the magnetic field, and the length, respectively, 
where is the density at the infinity, Hq = H z (0, 0) 
is the magnetic field at the center of our smoothed dipole 
field (equation 16), and R d is the radius of the current disk 
(equation 12). For the simulations presented here, the ra- 
tio of R d to the Bondi radius Rb = GM/c^ was chosen 
to be 



Rb 



1 



50V2 



(24) 



Note that for 7 = 7/5, the sonic radius of the Bondi flow 
is at r s = [(5 - 3 7 )/4]i? B = 10V2R d . 
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The values 



V A0 



Po 



Hi 
8tt 



(25) 



provide units for the velocity and pressure. After reduc- 
ing the equations (5)— (11) to dimensionless form, a non- 
standard plasma parameter, 







Ho 2 



24 
1 Vj 



(26) 



AO 



appears. This is the ratio of the plasma pressure at infinity 
to the magnetic pressure near the origin. This parameter 
connects the two parts of the considered problem, the gas 
dynamical unperturbed Bondi flow at large distances with 
the stellar dipole magnetic field at small distances. 

The Alfven radius Ra for quasi-spherical accretion onto 
a magnetic dipole is given roughly by the balance of 
ram pressure of the flow pv'jj <~ Mvff/(4wR 2 ) with the 
magnetic pressure H 2 /(87r) ~ p? /(8nR 6 ), where Vff = 
y/2GM/R is the free-fall speed. Thus we define 



R%= ( 



\2MV2GM 



2/7 



(27) 



Taking into account that for maximum Bondi accretion 
rate M = Mb = 47rA c -R 2 3 /9ooC o, we have 



r/th 



2/7 



where k = 



{R d /R B f /7 
(49V2 7 A C ) 2 / 7 



(28) 



Thus, our parameter from equation (26) determines the 
ratio of the Alfven radius (27) to the current disk ra- 
dius. The dependence R 1 ^ oc /3~ 2 / 7 for Rj/Rb fixed shows 
that j3 oc Mb I li ■ The important quantity Mb/ p 2 is re- 
ferred to as the "gravimagnetic" parameter by Davies & 
Pringle (1981) (see also Lipunov 1992). Equation (28) for 
the Alfven radius as a function of (3 is discussed further in 
§ 3.3. 

Another important dimensionless parameter of the 
model is the magnetic Reynolds number, 



Re, 



RpVo = WgpVb 

Vm C 2 



(29) 



where r\ m is the magnetic diffusivity, i?o is a character- 
istic scale and Vo is a characteristic speed of the accre- 
tion flow. The value of the magnetic Reynolds number 
Re m depends on the chosen length scale R . Here, it is 
appropriate to evaluate Re m at the Alfven radius where 
RqVq w \J2GMRa- Thus it is useful to introduce the 
dimensionless magnetic diffusivity, 



Vn 



1 



RaVa Re r 



(30) 



Numerical values are discussed in § 3.4. 

We do not offer a detailed explanation for the mag- 
netic diffusivity r\ m . However it could arise from three- 
dimensional MHD instabilities not accounted for in our 
two-dimensional (axisymmetric) simulations. An impor- 
tant instability is the interchange or Rayleigh-Taylor in- 
stability with wavevector in the azimuthal direction k = 



k^cj) (Arons & Lea 1976; Eisner & Lamb 1977). This in- 
stability allows blobs or filaments of plasma to fall inward 
across the magnetic field. Of course, perturbations with 
k<f, ^ are not allowed in the present axi-symmetric sim- 
ulations. 

2.5. Numerical Method 

Finite difference methods are used to solve the axisym- 
metric MHD equations (5) - (11). The calculations are 
done in five main stages: 

In the first stage, equations (9) and (10) for the fields 
A$ and H$ are solved with the right-hand sides of the 
equations set equal to zero. That is, we first solve for 
the advection of these fields. For this stage, we use the 
hybrid numerical scheme proposed by Kamenetskyi & Se- 
myonov (1989). This is based on two well known methods, 
the Lax - Wendroff method for the region where the solu- 
tion is "smooth" , and a numerical scheme with differences 
against the flow for the "non-smooth" regions. 

For the second stage, the finite conductivity is taken into 
account in the equations for and H^ with the advection 
terms omitted. That is, we now solve parabolic equations 
for these fields. For this, an explicit multistage numer- 
ical scheme is used. This scheme is called the "Method 
of Local Iterations" (MLI), and it is described in detail 
by Zhukov, Zabrodin, & Feodoritova (1993). The MLI 
scheme is explicit and absolutely stable. 

For the third stage, we solve the equations (5) - (11) 
including the magnetic terms but omitting the advection 
terms. At this stage the magnetic terms are known. 

For the fourth stage the full advection equations for p, 
pv, and pe are solved. For this we use an adaptation of 
Flux Corrected Transport method (FCT), based on ideas 
discussed by Boris & Book (1973). The transport of each 
density, for example, p, to the next time level is realized in 
several stages. In one of the stages explicit correction of 
the fluxes is done. The procedure is performed separately 
along the r and z coordinates using a dimensional split 
technique (Strang 1968). 

In the fifth and final stage, the Joule heating is taken 
into account in equation (11) for the internal energy. Over- 
all, our finite difference method is second order accurate 
in space and time for smooth flows. This method and nu- 
merical scheme was tested widely and also was successfully 
used earlier (Savelyev and Chechetkin 1994; Savelyev et al. 
1996, Toropin et al. 1997). 

3. RESULTS 

This work investigates spherical Bondi type accretion 
of plasma to a star with an aligned dipole magnetic field. 
The resistive MHD equations (5) - (11) were solved us- 
ing the method described in §2.5 and initial and boundary 
conditions described in §2.2 and §2.3. 

We performed simulations for 12 different values of 
(3 cx Mb/h 2 (equation 26), and magnetic diffusivities f] m 
(equation 30). The calculated flows in all runs are similar 
to that shown in Figure 2. In § 3.1, we describe in detail 
the nature of the flow for a representative run. In § 3.2 
we analyze the stationarity of the spherical accretion to 
a dipole. In § 3.3, we show the dependence of the flows 
on (3 and fj m . In § 3.4, sample runs of accretion to a ro- 



6 



TOROPIN ET AL. 




, V V v . N *X V 

, , , , \\ N N 
t \ \ \ \\ \\ \ 



-10 v 



////ftt\\\\\x\\ 



■ -2.5 



-10 



10 



■ 1.5 



log w (P) 



Fig. 2. — The figure shows an example of our calculated flow for accretion onto a non-rotating star with a dipole 
magnetic field at t — 2.5i //. The run is characterized by the dimensionless parameters (3 = 3.5 x 10~ 7 cx M^/fi 2 and 
magnetic diffusivity fj m — 10~ 5 , where Mb is the accretion rate and fi is the star's magnetic moment. The background 
gray scale represents the density of the flow and the solid lines the poloidal magnetic field lines. The length of the arrows 
(shown at every 32 th cell along both directions) is proportional to flow speed. The internal circular region represents 
the "accretor." The shock wave and associated transition from supersonic to subsonic flow is evident. The flow becomes 
strongly anisotropic close to the dipole. The simulations were done on a grid of 257 x 257 cells in one quadrant of the 
physical space. Only one quadrant is needed due to the assumed axisymmetry and the symmetry about the equatorial 
plane. 



tating star are presented, 
application of our results. 



In § 3.5, we give a numerical 



3.1. Illustrative Simulation Run 

Here we describe in detail a run with (3 = 3.5 x 10" 7 and 
fjm — Simulations were performed on a uniform grid 

with 257 x 257 square cells. The size of the computational 
region was 10Rd x WRd, the accretor radius was 0.5Rd, 
where Rd is the radius of the current disk. We measure 



time in units of the free 

,3/2 



-fall time from 



R„ 



= 0; 



that is, tff = RHL/V2GM. 

Spherical accretion to a magnetic dipole is very different 
from that to a non-magnetized star. Instead of supersonic 
steady inflow, which is observed in standard Bondi accre- 
tion, a shock wave forms around the dipole. The super- 
sonic inflow outside the shock becomes subsonic inside it. 
In all cases we observe that the shock wave gradually ex- 
pands outwards. Figure || shows the main features of the 




Fig. 3. — Enlarged view of the inner region of the plot of Figure ^. The anisotropy of the flow is evident. The torus-like 
region with small flow velocities is the "stagnation zone" . In the equatorial plane, the outer boundary of this zone is the 
magnetopause with radius R mp ~ 2.6Rd- There are two (± z) polar accretion columns. The outer dashed line represents 
the Alfven surface. At very small distance from the star, the flow in the polar columns becomes supersonic. The sonic 
surface is marked by a dashed line in the region of polar columns. 



SPHERICAL BONDI ACCRETION ONTO A MAGNETIC DIPOLE 



7 




Fig. 4. — The figure shows the radial variation of density p, radial velocity v r , temperature T, and ratio of total kinetic 
energy-density of the matter W — p(e + v 2 /2) to the magnetic energy-density E m = H 2 /87r in the equatorial plane 
(left column) for the acccretion flow presented at Figure || (t = 3.5i//). The positions of the Alfven radius r — Ra, the 
magnetopause r = R mp and the shock wave r = R s h are shown. Right columns show the z— dependences of the same 
parameters along the z— axis. The positions of the shock wave z = R s h and Alfven radius z = Ra are shown. 



10 20 30 40 50 60 70 SO 90 

6 (degrees) 



Fig. 5. — Differential mass accretion rate per unit solid angle dM/dQ as a function of the angle 9 (with respect to the 
z— axis) for spheres with radii Rd, 3Rd, 5Rd, and 9^ for the flow presented in Figure 2. At large distances (R = 5Rd) 
accretion is almost spherically symmetric. Closer to the dipole (R — 3Rd), the flow becomes anisotropic. At very small 
distances (R = Rd) most of the matter accretes to the poles along a narrow cone of half-angle 9 1 / 2 ~ 30°. 




Fig. 6. — Radial variation of the vertical magnetic field H z in the equatorial plane. The thin line shows the dependence 
of H z for a vacuum dipole field. The vertical dashed lines indicate the positions of the shock wave R s h, Alfven radius Ra, 
and magnetopause radius i? mp . 



flow at time t w 2.5£// when the shock has moved to the 
distance R s h — 8.lRd- We observe that for R > R s h the 
flow is unperturbed Bondi flow, whereas inside the shock 
for R < Rsh it is subsonic. Initially, the subsonic accre- 
tion to dipole is spherically symmetric, but closer to the 



dipole it becomes strongly anisotropic. Near the dipole 
matter moves along the magnetic field lines and accretes 
to the poles. Figure || shows the inner subsonic region 
of the flow in greater detail. The dashed line shows the 
Alfven surface, which we determine as the region where 
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Fig. 7.— Temporal evolution of the shock wave radius R s h, Alfvcn radius Ra, and magnetopause radius R mp - 
measured in units of tff which is the free-fall time from the distance R ma x- 
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Fig. 8.— (a) Mass accretion rate M through spheres of radii R at t\ — 2tff (dashed line) and ti = 2.8t// (solid line). The 
shock wave expands from position R = R s h(ti) to position Rshfo)- Inside the shock the flow is subsonic with accretion 
rate Mdip ~ 0.5Mb- (b) Time evolution of mass-fluxes M(R)/Mb throught spheres with radii R = Rd and R = 5Rd- 
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the matter energy-density W = p(e + v 2 /2) is equal to 
the magnetic energy-density E m = H 2 /(87r). The Alfven 
surface is ellipsoidal, with radius Ra = 3.1i?<j in the equa- 
torial plane, and Ra — 2.3Rd along the z— axis. Note, 
that the "theoretically" estimated Alfven radius (28) is 
i?^ 1 ps 2.6Rd- A significant deviation from spherically 
symmetric flow is observed for R <J 2Ra, because mag- 
netic field starts to influence the flow before it reaches 
the Alfven surface. Matter in the equatorial plane moves 
across the magnetic field lines, decelerates and stops at 
a radius r = R mp ss 2.6i?d, which we term the "mag- 
netopause radius." There is a torus shaped region - a 
stagnation region - which is avoided by accretting mat- 
ter. The flow velocities in this region are negligible. The 
magnetopause region is located inside the Alfven surface 
(see Figure ||). The accretion flow along the z— axis is ac- 
celerated and becomes supersonic at z ~ l.&Rd (see inner 
dashed line at Figure ||) . 

Figure ^ show the radial and axial variation of differ- 
ent parameters. One can see, that the density p is larger 
in the subsonic region compared with the Bondi solution. 
The velocity decreases by about a factor of 4 across the 
shock wave, and the temperature increases. The density in 
the magnetopause region R < R mp is lower than outside, 
while the temperature is larger (see Figure [|, left column). 
The bottom panel of Figure 4 shows the variation of ratio 
W/E m . The radius where W/E m = 1 is Alfven radius. In 
the z-direction (see Figure ||, right column) , matter moves 
along the magnetic field lines with the result that the flow 
parameters change smoothly. 

Figure |^ shows the transition from spherically sym- 
metric accretion outside the magnetosphere to highly 
anisotropic accretion along the polar columns within the 
magnetosphere. The figure shows the matter flux accreting 
through the unit solid angle dM/dtt at different inclina- 
tions 8 of this solid angle relative to the ±z axis. One can 
see that at large distances R = 5Rd, the accretion is al- 
most spherically symmetric, whereas at smaller distances 
it becomes more and more anisotropic. 

Figure || shows that the initial vacuum dipole magnetic 
field (Figure |l|) is strongly compressed by the incoming 
Bondi flow. As a result of interaction with Bondi flow, the 
magnetic field has a dipole dependence only inside Alfven 
radius, r <^ Ra, and it decreases faster than the dipole 
field for Ra < r < R s h- This is a result of an induced az- 
imuthal shielding current in the dipole's magnetosphere. 
Thus distributed current has a sign opposite to that of 
the star's intrinsic azimuthal current. For r J> R s h, the 
field decreases dramatically with r due to the shielding 
current. Thus, the magnetic field at the outer boundary 
of the simulation region is negligible. 

Figure [7| shows that the shock wave initially forms close 
to the magnetosphere at R ~ Ra and then gradually ex- 
pands outward. However, the Alfven radius Ra and mag- 
netopause radius R mp become steady after t ;> i///2 and 
remain steady thereafter. This means, that the magneto- 
sphere of the star reaches equilibrium with the incoming 
matter rapidly and this equilibrium does not change as a 
result of outward movement of the shock wave. 



3.2. Stationarity of Accretion Flow to the Dipole 



It is important to know whether the calculated accre- 
tion flow to the dipole is stationary or not. We analyze 
this here. 

After the passage of the shock wave, the new subsonic 
regime of accretion forms around the dipole. The region of 
subsonic flow expands together with the expanding shock 
wave. Here, we analyze the stationarity of the flow in 
the subsonic region. Figure ^a shows the distribution of 
matter fluxes through spheres of different radii M(R) as a 
function of R at t — 2tff and t = 2.8i//. One can see that 
matter flux outside the shock wave R > R s h is the Bondi 
rate Mb, while inside the shock wave it is significantly 
less. It is almost constant and equal to the accretion rate 
to the dipole M(R) m M dlp m 0.5M B - 

Figure |p shows the matter fluxes through fixed spheres 
located at radii R — Rd and R = 5Rd as a function of time. 
The matter flux through the sphere of radius R — Rd cor- 
responds to the matter flux to the dipole. It decreases and 
goes to the constant M = Mdi P at t > 1.5i//. The matter 
flux through the sphere R = 5Rd is initially the Bondi 
rate, but after passage of the shock wave it decreases to 
M = Mdi P and does not change thereafter. Thus, Fig- 
ures ||a and ||d demonstrate that the flow in the subsonic 
region is stationary in both space and time. The shock 
wave switches the Bondi flow to a new flow with stationary 
subsonic accretion and smaller accretion rate. The local 
physical variables, for example, density p and velocity v 
are also time independent. 

The formation of an expanding shock wave during accre- 
tion to a dipole results from the fact that the gravitating 
center with dipole field "absorbs" matter at a slower rate 
than the Bondi rate. The rate of accretion Mdi P at given 
parameters of the Bondi flow is determined by the physical 
parameters of the dipole (see § 3.3). 

An analogous situation was found in investigations of 
hydrodynamical accretion to a gravitating center. Kazh- 
dan & Lutskii (1977) (see also Sakashita 1974; Sakashita 
& Yokosawa 1974; Kazhdan & Murzina 1994) investigated 
spherically symmetric accretion flows for conditions where 
the matter flux through the inner boundary (which is the 
surface of the star) is less than matter flux supplied at the 
outer boundary. They found a family of self-similar so- 
lutions where the expanding shock wave links the regions 
inside and outside the shock wave. These regions have 
different stationary matter fluxes corresponding to matter 
fluxes at the boundaries. Our simulations show similar 
behaviour. 

Here, we should point out that the shock wave in our 
simulations is a temporary phenomenon, which establishes 
a new regime of accretion around the dipole. It appeared 
because the external accretion rate is larger than accre- 
tion rate which dipole can "absorb." A different situa- 
tion was considered by Ruffert (1994) who performed 3D 
hydrodynamical simulations of Bondi accretion. He used 
initial conditions where the matter distribution had con- 
stant density and zero velocity. In his simulations, the 
initial matter flux is zero and stationary accretion was es- 
tablished by a rareaction wave propagating outward from 
the central object. 

For fixed boundary conditions (supersonic Bondi ac- 
cretion at the outer boundary) the results of simulations 
do not depend on initial conditions. The dependence on 
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boundary conditions will be investigated separately. 

It is of interest to know how far outward the the shock 
wave will propagate. Note that the Bondi flow is super- 
sonic out to some distance and is subsonic at larger dis- 
tances. We expect that after reaching the subsonic area, 
the shock wave will "dissolve" and the flow will be purely 
subsonic. Thus, the shock wave may be only a temporary 
phenomenon which results from the initial conditions of 
our simulations. From the other side, if the flow is su- 
personic up to very large distances, then the shock wave 
expansion may be stopped by the physical structure of 
the astrophysical system. For example, in the wind-fed 
pulsars, the shock movement would stop at a radius of the 
order of the binary separation. 



3.3. Dependence of the Accretion Flow on (3 cx Mb/// 

and fj m 

We first analyze the dependence of the flow on the ex- 
ternal accretion rate Mb and the star's magnetic moment 
fj,. As discussed in § 2.4, these quantities are coupled so 
that the investigated physical model depends only on the 
ratio (3 cx Ms/// 2 . Each simulation run takes considerable 
time and for this reason we adopted the following proce- 
dure for deriving the dependence on (3. We start from the 
conditions of the simulation run of with = 10" 6 at time 
t = 2tff, and then change (3 by a factor of 10™/ 10 in a se- 
quence of 5 independent simulations (n = 1 ... 5). These 
simulations were performed up to t — 5i//. Fluctuations 
connected with the readjustment of the flow are damped 
by this time. We then measured the radius of the magne- 
topause R m p and the Alfven radius Ra for new values of 
(3. The Alfven radius in the equatorial plane is found to 
have a power law dependence, 



R A {(3) ra 2.3R d ^j^j , (31) 

where k ra 0.32 ± 0.03 (~2/Tw 0.286). The equatorial 
magnetopause radius is found to be proportional to the 
Alfven radius, R mp {(3) « 0.026i? d + 0.756R A {fl). Figure | 
shows the observed dependences. The Alfven radius is 
found to have a weak dependence on magnetic diffusivity 

— k 

rj m , Ra ~ ?7m , where k v ra 0.075. 

Figure O shows that the stationary accretion rate to 
the dipole Mdip also depends on (3. We find that M<u p is 
always smaller than the Bondi accretion rate Mb- The 



dependence found is 

M dip ra 0.78Mb (j^=e) " > ( 32 ) 

for (3 < 10~ 6 , where mp ra 0.51. 

Here, we recall that (3 cx Mb//* 2 and in equation (23) 
for Mb all quantities except are fixed in our simula- 
tions. Thus, the accretion rate to the dipole depends on 

3/2 

the density of surrounding matter as Mdip cx poo and on 
the star's magnetic moment as Mdip oc The last de- 

pendence can be explained by the fact that at larger (3, the 
Alfven radius is smaller so that the opening angle of the 
accretion columns is larger. As a result, matter can more 
readily flow into the gravitating center. At j3 > 3 • 10~ 6 , 
the dependence is different. Mdip increases more gradu- 
ally and approaches the critical Bondi accretion rate Mb- 
Here, we observe stationary Bondi flow similar to that ob- 
served in simulations to non-magnetized gravitating ob- 
ject (Ruffert 1994). 

Figure [l(] shows the dependence of the accretion rate on 
the magnetic diffusivity 

M dip ra 0.78M S (Jjpg) " , (33) 

for rj m < 10~ 5 , m v ra 0.38. This dependence means that 
matter accretes more readily at larger diffusivity as ex- 
pected. 

The half-opening angle #x/2 °f the accretion funnel at 
r = Rd (see Figure ||) decreases as fj m decreases. We find 

01/2 OC (fj m ) - 26 . 

At higher difusivities r) m > 10~ 5 , the dependence (33) 
becomes smoother, Mdip — ► Mb and we observe steady 
accretion at the Bondi rate. 

3.4. Accretion to a Rotating Dipole 

We also investigated cases of accretion to a rotating star 
with an aligned dipole magnetic field. To create the rotat- 
ing dipole, we rigidly rotated the current disk < r < Rd 
which is a part of boundary condition at z — as dis- 
cussed in §2.3. The disk radius R d is in effect the radius 
of the star. We discuss two cases, one with slow rota- 
tion, £1* = O.lf^if and the other with fast rotation, 
£1* = 0.35£l/< {Rd), where 11* is the star's angular rotation 
rate and Qj({Rd) — \JGM/R^ is the Keplerian angular 
velocity at the edge of the disk r — Rd- 




4 6 8 10 



P/IO" 7 

Fig. 9. — Dependence of the Alfven radius R A and magnetopause radius Rmp on parameter (3 cx Mb/ li 2 for magnetic 
diffusivity f) m = 10~ 5 . 
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Fig. 10. — Dependence of accretion rate to the dipole Mdi P , measured in Bondi accretion rate Mb, on j3 (at fj m = 10 5 , 
left frame) and on fj m (at j3 = 10~ 6 , right frame). 



We observed that in the case of slow rotation the gen- 
eral behavior of accretion flow is similar to that for the 
non-rotating case. The corotation radius R cr , where 
^k{Rct) = ^* or R cr = (GM/Ql) 1 / 3 , is significantly 
larger than the Alfven radius Ra for a non-rotating dipole. 

As in the non-rotating case, the shock wave forms and 
propagates outwards, while accretion to the dipole is sub- 
sonic and steady. However, a new feature appears: The 
stagnation region mentioned earlier rotates rigidly with 
the angular velocity of the star. We find that the limit 
of slow rotation is valid for R cr > Ra- In this case the 
linear velocity of rotation at the outer edge of magneto- 
sphere (R = Ra) is smaller than the Keplerian velocity. 
In cases of slow rotation, accretion to the dipole is steady 
but with an accretion rate M£°£ which is smaller than the 

corresponding value for a non-rotating star, Mdi P . 

The second case we discuss is a rapidly rotating dipole, 
where the corotation radius is smaller than the Alfven 
radius for the corresponding system with a non-rotating 
star, R cr < Ra- In this case the outer equatorial region 
of the rotating magnetosphere, R cr < r < Ra, has az- 
imuthal velocities in excess of Keplerian velocity. Mat- 
ter moves outwards in a wide "belt"- like region around 
the equatorial plane. Figures fTD and n2 show the simula- 
tion results for fl+ — 0.35f2/f when the corotation radius 
is R cr = 2.0i?d, and the Alfven radius in the equatorial 
plane is R r A ot = l.%Rd- One can see, that magnetic field 
lines are elongated in equatorial direction by outflowing 
matter. The strongest outward acceleration is in the equa- 
torial plane where the centrifugal force is largest. However, 
in the region of Alfven surface (see Figure P), an essen- 
tial acceleration is observed along the magnetic field lines. 
This acceleration appears to determine the unusual shape 
of the Alfven surface (see Figure [l2|). Also, an essential 
outflow may occur above and below the equator. 

At larger radial distances, the outflowing matter encoun- 
ters the incoming Bondi flow and turns to the direction of 
poles. The magnetic field lines are elongated in the r- 
direction. The Alfven radius in the direction of the poles 
has a value similar to that in the non-rotating case. How- 
ever, in the equatorial plane the Alfven radius is signifi- 
cantly smaller (i?™* « 1.9Rd) than its non-rotating value 
(Ra ~ 3-lRd)- Figure |l3| shows that the magnetosphere 
inside r < R r A ot rotates rigidly. The angular velocity de- 
creases gradually for r > R r A ot . The Figure also shows 
that the angular velocity of matter is larger than Keple- 
rian in the outer parts of the magnetosphere beyond the 
Alfven radius. In this region, the centrifugal force expells 
matter outward forming "propeller" - like outflows. The 
"propeller" outflows are predicted to occur in the mag- 



netosphere of a rapidly rotating magnetized star as first 
discussed by Illarionov and Sunyaev (1975). The theory 
of such outflows has received renewed interest for the case 
of disk accretion to rotating magnetized stars (Li & Wick- 
ramasinghe 1997; Lovelace, Romanova, and Bisnovatyi- 
Kogan 1998). A systematic study of spherical accretion to 
a rotating star with an aligned dipole magnetic field is in 
preparation (Toropin et al. 1998). 

3.5. Astrophysical Example 

Here, we present an application of our simulation results 
in terms of the physical quantities. We consider Bondi ac- 
cretion to a non-rotating magnetized protostar with mass 
M = 1 Mq = 2 x 10 33 g and radius R = 10 11 cm. We 
use our simulation run with f3 = 10 -6 and f) m — 10~ 5 , 
which is close to the case discussed in §3.1. We take the 
radius of the current disk to be equal to the radius of 
the star, Rd = R- From equation (24), the Bondi ra- 
dius is R B = 50V2R « 7.1 x 10 12 cm. The sound 
speed of the matter at infinity from equation (26) is 
Cqo = \/GM/Rb ~ 4.3 x 10 6 cm/s. This corresponds 
to a temperature Too w8x 10 4 K for a hydrogen plasma. 

We assume that the magnetic field at the star is Ho = 
100 G. Then, according to equation (16), magnetic mo- 
ment of the star is [i = 1.4 x 10 34 G cm 3 . From equations 
(26), the matter density at infinity is 

Poo = f£ ■ (34) 

Coo on 

For these parameters, the density at infinity is poo ~ 
3 x 10~ 17 g/cm 3 or a particle number density = 
1.8 x 10 7 l/cm 3 . The Bondi accretion rate from equation 
(23) is M B = 8.3 x 10~ 10 M Q /yr. 

Our code has finite magnetic diffusivity rj m . Here, we es- 
timate the magnetic Reynolds number Re rn . At a distance 
R from the origin we have 

Re m = ?^1L, (35) 



where vtf = \j2GMfR is the free-fall speed. Using the 
definition of dimensionless magnetic diffusivity, we get 
Vm = ^uRbVao = y / 2GMR B /j/3. Finally, we obtain 

« e ,„ = Mi(«)' = 

Vm \RbJ 
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Fig. 11. — The figure shows the accretion flow to a rapidly rotating star (see §3.3) at time t = 2.5tff . The gray-scale 
background represents the density and the solid lines the poloidal magnetic field. The length of arrows is proportional to 
the flow speed. 




Fig. 12. — Enlarged view of the accretion flow to a rapidly rotating star with dipole field (Figure |Tl|). The outer dashed 
line represents the Alfven surface. A inner sonic surface is indicated by the dashed line in the region of polar columns. 
The outflow in the equatorial plane starts from the region of the corotation radius R cr — 2Rd inside magnetosphere. 
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Fig. 13. — The radial variation of the angular velocity of the matter f2(r) = v,j,{r)/r in the equatorial plane measured 
in units of angular velocity of the disk (solid line). The thin solid line shows the radial dependence of the Keplerian 
angular velocity. The positions of the corotation radius R cr and R A ot are shown. Also, the position of the Alfven radius 
for the corresponding non-rotating system Ra is marked. 



At the distance of magnetopause, R — R mp ~ 3Rd, we get 
Re m ~ 2.9. 

4. CONCLUSIONS 

We have developed a method for MHD simulation of 
spherical Bondi-type accretion flow to a rotating star with 
an aligned dipole magnetic field. Using this method we 
have made a detailed study of the accretion to a non- 
rotating star for different accretion rates, stellar magnetic 
moments, and magnetic diffusivities. We also include an 
illustrative case of accretion to a rapidly rotating star. The 
simulation study confirms some of the predictions of the 
analytical models (Davidson & Ostriker 1973; Lamb et al. 
1973; Arons & Lea 1976). However, the simulated flows 
show a different behavior from the models in important 
respect summarized below. 

Our results for accretion to a non-rotating star agree 
qualitatively with some of the early theoretical predic- 
tions. In particular, (1) A shock wave forms around the 
dipole which acts as an obstacle for the accreting matter; 
(2) A closed inner magnetosphere forms where the mag- 
netic energy-density is larger than the matter energy den- 
sity; (3) The outer dipole magnetic field is strongly com- 
pressed by the incoming matter. (4) The flow is spherically 
symmetric at large distances, but becomes anisotropic near 
and within the Alfven surface. Closer to the star the accre- 
tion flow becomes highly anisotropic. Matter moves along 
the polar magnetic field lines forming funnel flows (David- 
son & Ostriker 1973). (5) The Alfven radius varies with 
[3 oc Mb/ as Ra ~ P~ a ' 3 , which is close to theoretical 
prediction Ra oc [3~ 2 ^ 7 (Davidson & Ostriker 1973). 

The new features observed in our simulations of accre- 
tion to a non-rotating star include the following: 

(1) We observe that the shock wave which initially forms 
around the magnetosphere is not stationary but rather ex- 



pands outward in all of our simulation runs. This is differ- 
ent from the theoretical models which assume a stationary 
or standing shock wave (Arons & Lea 1976); (2) A new sta- 
tionary regime of subsonic accretion forms around the star 
with dipole magnetic field; (3) A star accretes matter only 
at specific Mdi P rate which is less than the Bondi rate Mb- 
That is, Mdip = k Mb with k < 1; (4) This accretion rate 
Mdip is smaller when [3 oc Mb //J? is smaller, that is, when 
the star's magnetic field is larger. Also, Mdi P increases as 
the magnetic diffusivity r\ m increases. 

We are presently making a systematic study of accre- 
tion to a rotating star with dipole field. In this work 
we give only sample results which illustrate the new be- 
havior resulting from the star's rotation. Accretion to a 
slowly rotating star, where the corotation radius R cr = 
(Gm/Ql) 1 / 3, is significantly larger than the Alfven radius 
Ra is similar to accretion to a non-rotating star. How- 
ever, the rate of accretion Mdi P is smaller than in the 
corresponding non-rotating case. For a rapidly rotating 
star, where R cr < Ra, "propeller" outflows form in the 
outer parts of magnetosphere and outside magnetosphere 
as proposed by Illarionov and Sunyaev (1975). These out- 
flows result in a major change in accretion flow and field 
configuration. 
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